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ABSTRACT 


— J  The  general  planar  translation  of  two  bodies  of  revolution  through  an  inviscid  and 
incompressible  fluid  is  considered.  The  moving  trajectories  and  the  hydrodynamic  interactions  are 
computed  based  on  the  generalized  Lagrange's  equations  of  motion,  including  the  effects  of  solid 
constraints,  external  forces  in  the  plane  of  motion,  and  a  uniform  stream  in  any  direction  parallel  to  the 
plane  of  motion.  In  a  relative  coordinate  system  moving  with  the  stream,  the  kinetic  energy  of  the  fluid 
is  expressed  as  a  function  of  six  added  masses  due  to  motions  parallel  and  perpendicular  to  the  line 
joining  the  centers  of  the  solid  pair.  The  exact  solution  of  added  masses  in  closed  forms  are  obtained 
for  the  motion  of  two  spheres.  A  new  iterative  formula  based  on  the  analysis  of  velocity  potentials 
around  each  body  is  developed  for  added  masses  and  their  derivatives  with  respect  to  the  separation 
distance  due  to  the  transversal  motion.  The  method  of  successive  images  and  the  Taylor's  added-mass 
formula  are  applied  to  determine  the  added  masses  and  their  derivatives  due  to  the  centroidal  motion. 
These  results  are  compared  with  the  numerical  solution  of  added  masses  computed  by  the  boundary- 
integral  method  and  the  generalized  Taylor's  added-mass  formula.  The  integral  equations,  in  terms  of 
surface  source  distributions  on  both  surfaces,  are  carefully  modified  for  obtaining  accurate  numerical 
solutions.  Numerical  results  are  given  for  several  practical  engineering  problems.  ^  -  , 
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I.  INTRODUCTION 


The  hydrodynamic  interactions  among  solids  affect  the  motion  of  each  and  every  solid 
significantly  in  the  near  field.  If  the  Reynolds  number  based  on  the  size  and  velocity  of  a  typical  solid 
is  sufficiently  large,  the  inviscid  irrotational-flow  theory,  or  the  potential-flow  theory,  can  predict  the 
real  flow  with  sufficient  accuracy  and  the  effects  due  to  flow  separation,  boundary-layer  and  wake 
generation,  and  vortex  formation  may  be  neglected.  The  validity  of  this  assumption  has  been  shown 
by  Wu  and  Landweber*1*  by  comparing  the  prediction  with  the  measured  results  on  added  masses. 
The  potential-flow  theory  has  applications  in  a  variety  of  engineering  situations  such  as  the  impact  of 
floating  bodies,  the  motion  of  a  blunt  solid  around  other  fixed  or  moving  boundaries,  the 
hydrodynamic  interactions  on  bodies  due  to  the  oncoming  flow,  and  so  on. 

Hicks*2'  and  Herman*3*  first  analyzed  the  kinetic  energy  of  the  fluid  due  to  the  motion  of  two 
spheres  along  the  line  joining  their  centers,  and  obtained  analytic  solutions  of  added  masses  in  terms  of 
doublets  interior  to  each  body.  Their  expressions  about  the  strengths  and  positions  of  the  doublets 
were  alternatively  reduced  to  a  set  of  recurrence  formulas,  which  were  suitable  for  computation,  as 
shown  by  Landweber  in  the  book  edited  by  Rouse*4*.  The  evaluation  of  added  masses  due  to  the 
transverse  motion  is  more  complicated.  Hicks*2*  used  the  method  of  successive  images  again  and 
represented  the  added  masses  in  terms  of  distributed  and  isolated  dipoles;  however,  he  was  able  to 
calculate  only  a  few  images  owing  to  the  complexity  of  the  calculation.  Mitra*5*  and  Shail*6*  applied 
the  method  of  successive  images  to  the  calculation  of  potential  field  surrounding  both  spheres,  and 
obtained  an  analytic  solution  for  the  Dirichlet  problem  in  electrostatics.  They  took  two  sets  of 
spherical  polar  coordinates  at  the  centers  of  each  sphere  and  obtained  a  set  of  unknown  coefficients 
involved  in  the  series  expansion  of  velocity  potential  by  applying  the  Neumann-Liouville  iteration 
process.  By  changing  the  images  of  point  sources  to  dipoles,  we  can  extend  their  analysis  to 
determine  the  hydrodynamic  interactions  between  two  spheres  moving  along  their  centerline.  The 
numerical  solutions  of  the  potential  field  due  to  either  a  moving  spherical  body  in  the  inviscid  flow 
(Neumann  boundary-value  problem)  or  a  charged  spherical  body  (Dirichlet  boundary-value  problem) 
were  reported  in  detail  by  Atkinson*7*’*8*’*9*,  who  also  discussed  the  boundary-integral  approach  and 
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evaluated  various  solution  techniques  for  solving  the  integral  equations.  A  more  general  numerical 
model,  also  based  on  the  boundary-integral  method,  for  interaction  problems  of  two  bodies  has  been 
developed  by  Landweber  and  Chwang^10l  However,  their  results  indicate  that  the  numerical  error 
becomes  large  when  two  bodies  are  close  to  each  other.  The  motion  of  a  solid,  influenced  by 
hydrodynamic  interactions,  was  solved  by  Lamb^1 who  applied  the  Lagrange's  equations  of  motion 
in  the  generalized  coordinates  and  related  the  fluid  inertia  to  the  equations  of  motion  by  means  of  the 
kinetic  energy  of  the  fluid.  Recently,  Guo  and  Chwang^12^  applied  Lamb's  result  to  the  oblique 
motion  of  two  circular  cylinders.  They  also  modified  the  integral  equations  for  the  surface  source 
distributions  so  that,  when  two  circles  were  close  to  each  other,  the  steep  peak-values  in  the  integrands 
were  eliminated.  In  their  example,  very  good  agreement  of  the  numerical  result  with  the  exact  solution 
was  obtained  with  the  Gaussian  quadrature  formula. 

In  previous  studies  of  three-dimensional  solids  through  a  fluid,  much  attention  has  been 
focused  on  the  centroidal  motion  of  two  symmetric  bodies,  very  little  on  the  oblique  motion.  This  is 
mainly  due  to  difficulties  associated  with  the  added-mass  evaluation,  especially  when  two  bodies  are 
very  close  to  each  other.  Thus,  development  and  modifications  on  numerical  techniques  are  highly 
desirable  for  accurate  solutions. 

The  first  objective  of  the  present  work  is  to  consider  the  general  planar  translations  of  two 
bodies  of  revolution  which  are  symmetric  with  respect  to,  and  have  their  axes  of  rotation  perpendicular 
to,  the  plane  of  motion.  The  Lagrange's  equations  of  motion  will  be  generalized  for  the  trajectories  of 
moving  bodies,  including  the  effects  of  solid  constraints,  external  forces  in  the  plane  of  motion,  and  an 
unbounded  uniform  stream  in  any  direction  parallel  to  the  plane  of  motion.  The  second  objective  is  to 
obtain  the  exact  solution  of  added  masses  due  to  the  transverse  motion  of  two  spheres.  This  solution, 
together  with  that  for  the  centroidal  motion  of  two  spheres,  will  be  used  to  examine  the  reliability  of 
numerical  results.  The  third  objective  is  to  modify  the  integral  equations  which  govern  the  source 
distributions  on  each  solid  surface,  and  to  improve  the  accuracy  of  numerical  results.  We  shall 
consider  the  case  of  two  spheres  as  an  example  and  compare  the  numerical  result  with  the  exact 
solution. 
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The  generalization  of  equations  of  motion  for  two  bodies  in  translation  is  presented  in  Section 
2.  The  transformation  of  coordinates,  in  which  the  evaluation  of  added  masses  for  two  bodies  of 
revolution  may  become  simple,  is  also  given  in  this  section.  In  Section  3,  we  are  concerned  with 
analytical  solutions  of  added  masses  due  to  the  centroidal  and  transversal  motions  of  two  spheres.  We 
shall  develop  an  iterative  formula  to  calculate  the  unknown  coefficients  involved  in  the  series 
expansion  of  velocity  potentials  and  to  determine  the  related  added  masses.  The  numerical  solution  of 
added  masses  for  the  two-sphere  problem  is  given  in  Section  4.  The  set  of  integral  equations 
governing  the  source  distributions  on  the  surfaces  is  modified  in  order  to  improve  the  ill-behaved 
integration  when  two  bodies  are  close  to  each  other.  The  reliability  of  numerical  solutions  of  the 
integral  equations  obtained  by  the  Gauss-Seidel  iterative  method  and  the  Gaussian  quadrature  formula 
is  examined  by  comparing  them  with  the  analytic  solutions.  Several  examples  are  given  and  discussed 
in  Section  5.  Finally,  conclusions  are  presented  in  Section  6. 


II.  EQUATIONS  OF  MOTION 

Consider  two  solids  translating  in  an  unbounded,  inviscid,  and  incompressible  fluid  which 
moves  with  a  uniform  velocity  UQ  at  infinity.  Let  xa  and  Ua  denote  the  instantaneous  position  and 

velocity  vectors  of  body  a  (a  =  1,2)  in  a  relative  coordinate  system  moving  with  the  uniform  flow; 
that  is, 

Ua  =  u<x-Uo’  (D 

where  ua  is  the  absolute  velocity  of  body  a. 

From  Lamb^1  the  force  acting  on  body  a  by  the  fluid,  due  to  hydrodynamic  interactions,  is 
governed  by  Lagrange's  equations  of  motion  in  generalized  coordinates. 


F  -  d  aT  il_ 

* "'**»■.  a*' 


(2) 


where  t  is  the  time,  Uio  and  xia  denote  the  i-th  (i  =  1,  2,  3  in  general)  components  of  Ua  and  xa 
respectively,  and  T  is  the  kinetic  energy  of  the  fluid.  If  there  is  an  external  force  Eja  acting  on  body  a 
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(3) 


of  mass  Ma  in  addition  to  the  force  due  to  the  fluid  pressure,  the  equations  of  motion  become 


E  .  +  F  -  M 
^ia  +  ria  ~  Ma  dt 


(no  sum  on  a) 


The  kinetic  energy  T  in  equation  (2)  is  given  by  the  integration  over  the  solid  surfaces, 


2T  =  -p  ^dS  , 


where  <(>  is  the  velocity  potential,  n  is  the  unit  outward  normal  of  each  surface  and  p  is  the  fluid 
density.  In  the  relative  coordinates  moving  with  the  uniform  stream  at  infinity,  the  fluid  is  at  rest  at 
infinity  and  disturbed  only  by  motions  of  bodies,  and  the  flow  is  therefore  irrotational.  The  velocity 
potential  <j>  satisfies  the  Laplace  equation  and  may  be  expressed  as 

♦  =  *ia  Uia>  (5) 

where  <|>ia  is  the  velocity  potential  due  to  the  unit  velocity  of  body  a  in  the  i-th  direction.  Summation 
on  the  repeated  indices  is  implied  unless  indicated  otherwise.  From  equations  (4)  and  (5),  the  kinetic 
energy  can  be  expressed  as  (see  Landweber  and  Chwang^10') 


2T  ■  A»jP  utaujP 


(i,  j  =  1 .2,3  and  a,  p  =  1,2), 


where 


Aiaj$  =  Ajpia  =  -P  4>i«  ^  dS  (no  sum  on  p) 


are  added  masses.  Referring  to  the  expression  of  hydrodynamic  forces  given  by  Landweber  et  al } 13', 
we  can  describe  the  general  translations  of  two  bodies  in  the  relative  coordinates  by 

MaTt“  =  Eic  *  Aiajp^  •  (Uki-Uid)  V  \  UjPU"*^mi  <5«1'  5«2> 


(no  sum  on  a)  ,  (8) 

where  sk  (k=l,2,3  in  general)  is  the  k-th  component  of  (xj  -  x2).  There  are  six  differential  equations 
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in  (8)  for  two  bodies  translating  in  three  dimensions.  The  initial  conditions  required  for  solutions  of 
(8)  are  given  based  on  physical  problems. 

The  first  integration  of  equation  (8)  yields  either  the  velocity  components  of  each  solid  relative 
to  the  uniform  flow  at  any  time  instant,  or  the  unknown  forces  acting  on  bodies  due  to  solid 
constraints;  and  the  second  integration  gives  the  trajectories  of  each  body.  Since  the  rotational  motion 
is  not  considered  in  the  present  derivation,  the  equations  of  motion  (8)  is  applicable  to  (i)  the 
translation  of  two  spheres  in  three  dimensions,  (ii)  the  planar  translation  of  two  bodies  of  revolution 
which  are  symmetric  with  respect  to,  and  have  their  rotating  axes  perpendicular  to,  the  plane  of 
motion,  and  (iii)  the  centroidal  translation  of  two  bodies  which  are  symmetric  with  respect  to  the 
centerline.  In  these  three  cases,  the  rotation  of  each  body  vanishes  since  the  moments  due  to  the 
hydrodynamic  interactions  are  zero.  However,  in  the  present  study,  we  will  focus  our  attention  on  the 
first  two  cases. 

For  the  planar  translation  of  two  bodies  in  the  xy-plane  (Fig.  1),  we  shall  use  the  simpler 
notation  (i,  j  =  1,2) 


uu=  uit 

Eu=Ei,  Ui2 

-  Ui+2,  Ei2  -  Ej+2, 

Ailjl  =  Aij’ 

Ailj2  =  Ai(j+2)’ 

Ai2j2  =  A(i+2)(j+2) 

Thus,  equation  (8)  reduces  to 


(9a) 

(9b) 


Mi^  =  Er^-U/Uk-Uk.2)^  +  J^(\rWUiUm  ("O  sum  on  i),  (10) 

where  k  =  1,2  and  other  subscripts  have  a  range  of  1  to  4,  Mi  =  M2  is  the  mass  of  body  1  and  M3  = 


dA:;  1  9  A;, 


M4  mass  of  body  2. 


The  added  masses  appearing  in  equation  (10)  are  evaluated  in  the  Cartesian  coordinates  shown 
in  Fig.  1 .  For  two  bodies  of  revolution  which  are  symmetric  with  respect  to,  and  have  their  rotating 
axes  perpendicular  to,  the  xy-plane.  the  evaluation  of  added  masses  may  be  simplified  significantly  in 
another  Cartesian  coordinates  (x',  y'),  where  the  x'  axis  is  parallel  to  the  line  joining  the  centroids 
(Fig.2).  In  this  new  coordinate  system,  the  translations  of  o,  and  o2  can  be  decomposed  into 
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components  along  the  x'  and  y'  axes  and  the  new  added  masses  A'y  (i,  j  =  1,  2,  3,  4),  which  are 

related  to  the  centroidal  or  transversal  motion  of  bobies,  become  functions  of  the  separation  distance  s 
only.  Since  body  1  and  body  2  are  symmetric  with  respect  to  the  centerline  o^,  a  sign  change  in  U'2 
should  not  change  the  fluid  kinetic  energy  in  terms  of  A’y.  Therefore  A'12  =  A'32  =  0.  Similarly, 
A’14  =  A'34  =  0.  Thus,  the  fluid  kinetic  energy  T  in  this  coordinate  system  is  reduced  to 


2T  =  A^U*  +  2A'13U\U3  +  A33U323  +  A^U*  +  2A^U2U4  +  A^ 


Let  y  be  the  angle  between  the  \l  axis  and  the  x'j  axis  (Fig.2).  The  transformation  of  velocities 
is  given  by 


U'i-byUj,  (12) 

where  by  is  the  transformation  tensor  and  its  matrix  representation  is 

cosy  siny  0  0 

-siny  cosy  0  0 

0  0  cosy  siny 

0  0  -siny  cosy 

Zeros  in  B  indicate  that  there  is  no  constraint  between  body  1  and  body  2.  Based  on  the  invariance  of 
the  kinetic  energy  to  the  coordinate  transformation,  we  have 


The  added  masses  A'y  are  functions  of  the  separation  distance  s  only,  and  by  are  functions  of  the  angle 
y.  Equation  (14)  can  also  be  expressed  in  matrix  form  as 

A=Bta'B,  (14b) 

where  A  =  [Ay]  and  A'  =  [A'y]  are  matrices  of  added  masses.  By  (13)  and  (14),  we  can  explicitly 
write  down  the  added  masses  Ay  in  terms  of  the  added  masses  A’y  due  to  the  centroidal  and 
transversal  motions  with  respect  to  the  centerline.  Thus 
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An  =  A'ncos2?  +  A'22  sin2? ,  A13  =  A'^cos2?  +  A^sin2? , 

A12  =  (A'n  -  A22)  sin?  cos? ,  Au  =  (A'n -  A'24)siny cosy , 

A22  =  A' 11  sin2?  +  A22  cos2?  ,  A23  =  (A' 13  -  A'24)sin?  cos? , 

A34  -  (A33  -  A’44)  sin?  cos? ,  A33  =  A33COS2?  +  A44  sin2? , 

A24  =  A'nsin2?  +  A'24  cos2? ,  A44  =  A^sin2?  +  A’44  cos2? ,  (15) 


The  derivatives  of  added  masses  with  respect  to  the  k-th  component  of  the  separation  distance  sk  (k=l, 


2)  are 


3A™  h  h  aAii  9s  9bim  3bin  9? 


9An  9An  sin?9An 

— u  =  cos?  — 11  -  — 1  — — 
3sj  9s  s  9? 


dA, ,  0A1 ,  cos?  oA  1 1 

— 11  =  sin?  1  — 11 ,  etc. 

9s2  9s  s  9? 


where  Sj  and  s2  are  the  components  of  (xj  -  x2)  in  the  x  and  y  directions  respectively. 


(16b) 


Equations  of  motion  (10)  have  to  be  decoupled  before  we  can  solve  them  numerically.  Let 


G(U)  =  (U,  ■  U3)  Ajj  +  (U2  -  U4)  A;2,  (17) 

where  ^  =  9A/9sk  (k=l,2)  are  matrices  formed  by  the  derivatives  of  added  masses  (16),  and  M  is  a 
4x4  diagonal  matrix  with  elements  Mj,  M2,  M3,  M4.  Let  q(U)  be  a  column  vector, 

q(U)  =  \  [UT  ASj  U,  UT  Ag2  U,  -UT  \  U,  -UT  U]T.  (18) 


Then  equation  (10)  may  be  written  in  vector  form  as 

^-  =  (M  +  A)-1  [E  -  G(U)U  +  q(U)],  (19) 

where  E  =  [Ej,  E2,  E3,  E4]  is  a  column  vector  representing  external  forces  acting  on  bodies  1  and  2. 
Replacing  U  by  (u  -  U0),  where  U0  =  [U0l,  Uo2,  Uol,  Uo2]T,  Uol  and  Uo2  are  velocity  components 
of  the  uniform  flow  in  the  x  and  y  directions  respectively,  we  obtain  from  (19)  the  equations  of  motion 
in  the  absolute  coordinates. 
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(20) 


^  =  (M  +  A)'1  [E  -  G(u)u  +  G(u)U0  +  q(u-U0)]. 

There  are  various  numerical  techniques  for  solving  the  initial-value  system  presented  by 
equation  (20).  In  the  present  study,  the  Runge-Kutta-Fehlberg  method,  which  was  discussed  in  detail 
by  Atkinson^14!,  are  used  for  the  solution  of  velocity  components.  The  size  of  the  time  step  is  adaptive 
based  on  a  pre-assigned  error-control  parameter  in  the  calculation. 

III.  EVALUATION  OF  ADDED  MASSES  AND  THEIR  DERIVATIVES 

If  the  velocity  potential  due  to  the  motion  of  solids  can  be  represented  by  a  set  of  isolated  or 
distributed  sigularities  interior  to  solids,  the  well  known  Taylor’s  added-mass  theorem  and  its 
generalization  were  recommended  by  Landweber  and  Yih^15^  to  determine  the  added  masses.  For  a 
pair  of  three-dimensional  solids  moving  in  any  manner  in  the  xy-plane  except  pure  rotations,  the 
Taylor’s  added-mass  formula  can  be  generalized  to  (see  Guo  and  Chwang^12') 

Ay  +  5ij  M'j  =  4np  [  j  -xj  X;dV  +  £  (n^  x°j  +  p.y)  ]  (no  sum  on  j),  (21 ) 

f 

where  Mj  is  the  mass  of  the  fluid  displaced  by  body  1  (for  j  =  1.2)  or  body  2  (for  j  =  3,4),  Xj  and  m; 

are  the  volume-distributed  source  density  and  the  isolated  source  strength,  respectively,  inside  body  j 
due  to  the  i-th  velocity  component,  p^  is  the  strength  of  an  isolated  dipole  in  the  j-th  direction  inside 
body  j  associated  with  the  i-th  velocit>  component,  Xj  and  x°j  are  the  j-th  local  coordinate  of  X;  and  irij 
respectively  with  respect  to  body  j,  the  integration  is  over  the  volume  of  body  j,  Vj,  and  the  summation 
is  over  all  isolated  singularities  inside  Vj. 

In  the  case  of  centroidal  motion  of  two  spheres  along  their  centerline,  it  is  well  known  that  the 
velocity  potential  of  the  fluid  due  to  the  unit  motion  of  each  sphere  can  be  simply  represented  by  a  set 
of  isolated  doublets  inside  spheres.  The  location  and  strength  of  each  doublet  are  determined  by  the 
sphere  theorem  which  states  that  an  isolated  doublet  of  strength  p.  at  a  point  P  outside  a  sphere  of 
radius  a,  pointing  along  the  radial  axis  of  the  sphere,  has  its  isolated-doublet  image  of  strength  p(aA)3 
at  the  inverse  point  inside  the  sphere  with  the  direction  opposite  to  the  original  one,  where  X  (X>a)  is 
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the  distance  between  the  point  P  and  the  center  of  the  sphere.  Guo  &  Chwang^12'  have  discussed  the 
strengths  and  locations  of  the  image  doublets  for  the  motion  of  two  cylinders  and  derived  expressions 
of  added  masses  and  their  derivatives  by  using  Taylor's  added-mass  formula.  With  slight 
modification,  their  result  can  be  extended  to  added  masses  and  their  derivatives  due  to  centroidal 
motion  of  two  spheres.  By  applying  the  Taylor's  added-mass  formula  and  the  sphere  theorem  to  the 


centroidal  motion  of  two  spheres,  we  have 


Au=4jtp[  I  -  y  1, 
n=0 


A  33  —  4ttp  [  X  "  3  ]» 
n=0 


a  i3 =  4ttP  m2„+i’ 
n=0 

A31  =  A13’ 


where  is  the  strength  of  the  n-th  image  doublet  inside  sphere  1  of  radius  a  due  to  the  unit  motion  of 
sphere  1  along  the  x'  direction  (U'j  =  1,  U'2  =  U'3  =  U'4  =  0),  |J.2n+1  *s  lhe  strength  of  the  image 
doublet  of  p^n  inside  sphere  2  of  radius  b,  p*2n  denotes  the  strength  of  the  n-th  image  doublet  inside 
sphere  2  due  to  the  unit  motion  of  sphere  2  in  the  x'  direction  (U'3  =  1,  U'j  =  U'2  =  U’4  =  0),  and  p 

is  the  density  of  the  fluid.  Analogous  to  the  analysis  of  Guo  &  Chwang1121,  we  can  show  that 
sequences  {p2n}  and  (M-2n+iJ  ^  uniformly  convergent  for  s  in  the  region  [a+b,°°),  while  their 
derivatives  with  respect  to  s  are  uniformly  convergent  in  the  region  (a+b,°°)  but  divergent  as  s 
approaches  (a+b).  When  s  >  (a+b),  the  derivatives  of  added  masses  are  given  by 


ds  P  \  ds 


^33  ;  d^2n 

ds  P  ds 


iA13  .  ^  d^2n+l 

ds  P  ds  ’ 


dA3]  dAj3 


The  general  strengths  and  locations  of  image  doublets  and  their  derivatives  with  respect  to  s  are 


K  =  s. 


— 2  =  l 
ds  l’ 


determined  by  the  iterative  formula, 

M»=a3«. 
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where  denotes  the  distance  form  o2  to  the  n-th  inverse  point  in  sphere  1  and  X2n+1that  from  Oj  to 
the  (n+l)-th  inverse  point  in  sphere  2.  All  inverse  points  lie  on  the  centerline  0j02.  The  value  of  (i  2n 
is  obtained  directly  by  interchanging  a  and  b  in  equation  (24).  The  limiting  values  of  added  masses, 
as  s  approaches  (a+b),  are  derived  as 


where  £(3)  is  the  zeta  function.  As  s  tends  to  (a+b),  however,  the  n-th  terms  of  sequences  (d|i2n/ds) 
and  {dM2n+1/ds}  do  not  approach  zero  as  n  goes  to  infinity,  and  the  derivatives  of  added  masses  are 
divergent.  Thus,  when  two  spheres  are  very  close  to  each  other,  the  kinetic  energy  of  the  fluid  due  to 
the  centroidal  motion  of  these  two  spheres  is  finite  but  the  hydrodynamic  interaction  forces  approach 
infinity. 

When  two  spheres  make  transversal  motion  perpendicular  to  the  centerline  o,o2,  the 

determination  of  the  strengths  and  locations  of  the  hydrodynamic  singularities  become  very  difficult. 
Consider  the  unit  motion  of  sphere  1,U2  =  1,U1=U3  =  U4  =  0.  If  sphere  2  were  absent,  the 
velocity  potential  due  to  this  motion  could  be  represented  by  an  isolated  doublet  located  at  Oj  in  the 
direction  of  U  2.  However,  the  presence  of  sphere  2  violates  the  boundary  condition  on  surface  2  and 
it  requires  images  of  the  isolated  doublet  to  satisfy  the  boundary  condition.  Since  this  doublet  is 
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perpendicular  to  the  radial  axis  of  sphere  2,  its  images  include  an  isolated  doublet  at  the  inverse  point 
inside  sphere  2  plus  a  line  distribution  of  doublets  from  the  inverse  point  to  the  center  of  sphere  2.  The 
isolated  and  the  line-distributed  doublets  in  sphere  2  have  another  set  of  images  in  sphere  1,  which 
include  an  isolated  doublet  at  the  inverse  point  and  a  more  complicated  line  distribution  of  doublets 
from  the  inverse  point  to  Oj.  Very  rapidly,  the  expression  of  image  doublets  in  each  sphere  becomes 

extremely  complicated.  Hicks ^  stopped  his  calculation  at  the  fourth  image  system  because  of  the 
exceedingly  laborious  work  and  presented  an  approximation  of  the  velocity  potential  due  to  the 
transversal  motion  when  the  separation  distance  s  is  considerably  large. 

In  the  present  study,  however,  we  shall  consider  directly  the  velocity  potential  of  the  fluid 
instead  of  the  hydrodynamic  singularities.  Let  us  consider  two  sets  of  spherical  coordinate  systems 
(r1,91,X1)  and  (r2,02,X^)  with  origins  located  at  Oj  and  o2  respectively,  the  common  polar  axis  being 
0j02  (see  Fig.  3).  Let  <|>(l)  be  the  velocity  potential  due  to  the  unit  motion  of  sphere  1  in  the  direction  of 

U'2  and  expressed  in  terms  of  the  i-th  (i  =  1,2)  set  of  the  spherical  coordinate  system.  In  the  absence 
of  sphere  2,  \  the  velocity  potential  due  to  the  motion  of  sphere  1,  would  be  the  solution  of  the 


Laplace  equation  satisfying  the  boundary  condition  on  surface  1.  If  sphere  2  is  inserted  into  the  field  at 
o2,  we  should  add  an  extra  term  to  the  velocity  potential  <t>Q2),  which  is  (j)^  expressed  in  the 


(r2,02,X.2)  coordinate  system,  in  order  to  satisfy  the  boundary  condition  on  surface  2.  Thus,  the  new 
expression  for  the  velocity  potential  around  these  two  spheres  becomes 


*<2)=  *?>+♦?’ 


n\ 

However,  the  added  term  <t>j  violates  the  boundary  condition  on  the  surface  1  again  and  requires  an 
additional  term  ^  to  satisfy  the  boundary  condition  on  surface  1.  Thus 


♦(I> = C + oi” + 4". 


Continuing  this  process,  we  can  determine  the  velocity  potential  due  to  the  transversal  motion  of 
sphere  1  in  the  y'  direction  (U'2  =1,  U'j=  U'3=  U’4=  0)  as 
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(26a) 


♦0)=  z  ♦!!’ . 


►®  -  r • 

m=0 


(26b) 


In  the  spherical  coordinates  (r,,  0j,  Xj),  <j)(1)  satisfies  the  Laplace  equation  and  the  boundary 


condition  on  sphere  1, 


(VO  0, 
where  X  =  Xj  =  X^  and 


^’’(a^A)  .  „ 

- - —  =  sin  0,  cos  A., 

8r, 


13,13,  1  3  .  .  .  3  ,  1  92 


(y(0)2  = - (r?  — )  + - (sin0| — )  + - 

.?*.  Ame.39,  '39,  J  sin2e. 

t  1  1  1  1 


5?  (i  =  U) 


are  the  Laplacian  operator  in  (r;,  8j,  X).  In  the  coordinates  (r2, 02,  X),  <J(2)  satisfies 


(V(2))2  <|><2>  =  0, 


a<t>(2)(b,02,X) 


Around  sphere  1,  the  general  solution  of  equation  (27)  is  given  by 
♦“>  =  I  - +  Z  An(s)  R„  (r,"+  STT  ~ |  )1  T  cos  X, 


ri  n=l 


and  in  the  neighborhood  of  sphere  2,  the  solution  of  equation  (28)  is 


=  ( y  cos  X)  ZBn(s)Sn(r2"+sL^), 


(29b) 


where  Rn  =  P*(cos  00  is  the  associated  Legendre  function  of  degree  n  and  of  order  1 ,  S„  =  P*  (cos 


02),  and  Bn  are  arbitary  coefficients  which  depend  on  the  separation  distance  s  and  should  be 
determined  by  the  boundary  conditions  on  both  surfaces. 
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The  transformation  of  and  Sn  in  these  two  sets  of  spherical  coordinate  systems  (r1§  0j,  X) 
and  (r2, 02,  X)  will  be  used  frequently  in  the  following  analysis.  For  the  associated  Legendre  function 
of  degree  n  and  of  order  m,  the  general  transformation  was  first  given  by  Basset^16^ 


P^cosej) 

-  n+1 


(-1  )n 


m 


(n-m)!  sn+m+1  k_Q  (2m+k)! 


_  (n+m+k)!/T k pm  \ 

2-  f2m+ki!V  Pm+k(cos02)’ 


(30a) 


pj1(cos02) 

,  n+1 
r2 


m 


(n-m)!  sn+m+1  k_Q  (2m+k)! 


v  (-l)k(n+m+k)!  Ak  m  ,  , 

£  - r2m+kl!  (s}  Pm+k(cos0i)- 


(30b) 


For  m  =  1,  the  transformation  is  simplified  to 


_ (lU 


n-1 


(n+k)!  r2  k 
r1n+1  ”  (n-1)!  sn+1  ^  <k+l)!  (s} 

Sn  1  ”  (-l)kjl(n+k)l  /!  xkn 

r-1=(n-l)!sn+1  Cl  <k+1>!  <s)  k‘ 


(31a) 


(31b) 


By  the  ratio  test,  we  obtain  that  for  any  fixed  n,  the  radius  of  convergence  of  the  summation  in  (31a)  is 
r2  <  s  and  that  of  the  summation  in  (3  lb)  is  r,  <  s. 

If  sphere  2  were  absent,  from  (30a),  the  velocity  potential  would  be 

.  ,3 

L(! 

rI 

From  the  transformation  (31a),  can  be  expressed  in  the  second  coordinate  system  as 


\>o]  =  -^  (y  cosX). 


<2)  =  (yCOsX.)  I  Aln1JySn  , 
n=l 


HI  n 


where  A^1*  =  -  s'^n+2\  By  comparing  the  velocity  potential  <j>g2)  with  the  general  solution  (29b),  we 
should  add  the  term 
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,3  00  n  b2n+1 

n+1  _  n+l 


<>i2)  =  (%-  cosX)  I  A[n1]Sn 


n=l 


to  4>q2)  in  order  to  satisfy  the  boundary  condition  on  surface  2.  Thus 


„<2)  A  a<2)  _  /* 


+  f"  =  (%cosX)  lAf's.d/tjJj 


n  b2n+1 


n=l 


ri' 


n+l 


), 


which  is  convergent  as  r2  <  s.  The  term  <|>j2),  which  can  be  expressed  in  the  first  coordinate  system  as 

4>(11)  =  (ycosX)  E  A#21  fj  Rj,, 
n=l 


where 


(-l)"'1  I 

k=l 


k2  b2k+1  (n+k)!  [i] 

(n+l)!  (k+1)!  sn+k+1  k  * 


makes  an  extra  contribution  to  the  velocity  potential  and  in  turn  violates  the  boundary  condition  on 
surface  1.  This  contribution  should  be  corrected  by  adding  ^  to  the  velocity  potential  <)>j  , 


4>(i1)  +  ^1)  =  (ycosX)  I  A(n21Rn(r" 
n=l 


n  a 


2n+l 


(n+l)  r 


n+l 
1 


0. 


which  satisfies  the  boundary  condition  on  sphere  1  but  not  on  sphere  2.  The  summation  is  convergent 
as  rt  <  s.  Repeating  the  same  process  indefinitely,  we  obtain  the  velocity  potential  <)>(1)  and  <J)(2\  given 

by  (29a)  and  (29b)  respectively,  with  the  coefficients 


An(s)  = 


IA 


[2m] 

n 


m=0 


Bn(s)  =  S  A 
m=0 


[2m+l] 


where 


A[l)_  _-(n+2) 


(-l)"'1  I 


k=l 


k2  b2k+1  (n+k)!  ^[2m-i] 

(n+l)!  (k+1)!  sn+k+1  k 


(32) 
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.  [2m+i]  _  y  (-l)k-1k2  a2k+1  (n+k)!  r2m] 
^  k=l  (n+1>!  (k+1)!  sn+k+1  ^ 


(33) 


It  can  be  shown  that  both  Ajj(s)  and  Bn(s)  are  less  than  Cs‘n,  where  C  is  a  positive  constant. 
Therefore,  <|>^  and  are  convergent  uniformly  in  the  neighborhood  of  sphere  1  and  sphere  2 
respectively. 

By  equations  (7),  (9),  (27),  and  (29),  and  the  relation  <j>(1)  =  we  have 
2n  k  oo  oo 

A*22  =-  Jcos2Xd\  Jsin20,[  -  ^  +  I  (  IA[2ral)  ^  anRn]d0i 

O  O  a  n=i  m=0 


:  -  M*i(l  -  -  X  Ai2m^a3) . 
2  2m=0 


(34a) 


In  the  integration,  all  terms  containing  Rn,  n>l,  vanish  because  of  the  orthogonality  relation. 
Similarly,  A  24,  and  A  M  are  given  by 

27C  k 


A 24=  -  Jsin20,cos2>.  I  (  SA^^Sjd©,*-^^2  I  A1,2"1"11, 

n_  j  m=0  m=0 


OO  OO 


00 


0  o 


(34b) 


A„  =  -^  JdX  Jsin202Cos2M  -7T  +  £  ( I  Af”1)  b"S„  ]d02 


2k  n 


Soo  OO  ~  - 

1  /  »-*  »  *[2m].  211+1  .  nr 


b2  n=1  m=0 


_  1  ly|'  /  ■«  3b  a  *[2m], 

2M2(I'Tj0A1  >' 


where  A*[2ml  are  calculated  by  interchanging  a  and  b  in  the  iterative  formula  (33). 


(34c) 


The  derivatives  of  added  masses  with  respect  to  the  separation  distance  s  are  simply  given  by 

A[2m+1) 

dAH  3M,a3  ”  ' 

ds  d  *■ 


4  nSo  ds  ’ 


dA,4  3  M'?a3  “ 


ds 


4  mto  ds  ’ 
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oo 


(35) 


dA 


ds 


44  _  3M2b3  “ 
c  *  4  ~ 


dA! 


[2m] 


m=0 


ds 


From  formula  (33),  the  derivatives  of  and  Aj[2m+1l  with  respect  to  s  are 


dA™ 

ds 


=  0, 


ds 


=  (n+2)  s‘(n+3). 


dA^2  dA^21  (n+k+l)A[211 

_ _ _ (  n„-i  t?  k2b2k+1(n+k)!  1  2 _ _ 

ds  L  (n+l)!(k+l)!  lsn+k+t  ds  *  sn+k+2 

iC~"  1 


], 


dAi2m+11  oo 

=  I 


(-lV^a^Qi+k)!  r  1 


dA[2m)  (n+k+1)  A*2"11 


ds 


k=l 


(n+l)!(k+l)!  ^n+k+t  ds 


gn+k+2 


(36) 


We  shall  define  the  added-mass  coefficients  by 


kjj  =  A  ij  /  M  2 


(37) 


For  the  centroidal  and  transversal  motions  of  two  equal  spheres,  a/b  =  1,  the  values  of  ky  are  given  in 
Table  1  with  the  separation  distance  s  varying  form  2.01  to  10. 


IV.  NUMERICAL  SOLUTION  OF  ADDED  MASSES 
In  order  to  extend  the  above-mentioned  analysis  for  a  pair  of  spheres  to  the  general  case  of  two 
bodies  of  revolution,  we  may  have  recourse  to  the  numerical  computation  of  added  masses,  since  the 
exact  formulas  such  as  equation  (34)  do  not  exist  in  general.  In  a  given  geometric  and  kinematic  state, 
the  Neumann  boundary-value  problem  of  two  bodies  needs  to  be  solved  numerically  for  either  the 
velocity  potential  or  the  strengths  of  the  surface  distribution  of  sources  on  each  body.  We  shall 
formulate  the  problem  by  the  well-known  boundary-integral  method,  which  leads  to  a  pair  of 
Fredholm  integral  equations  of  the  second  kind,  and  consider  the  reliability  of  numerical  solutions 
obtained  from  different  procedures. 
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IV.  1  Formulation  of  integral  equations 

The  boundary- integral  method,  derived  from  the  generalized  Taylor's  formula  and  the 
fundamental  relations  between  normal  velocities  and  velocity  potentials  on  solid  surfaces,  has  been 
studied  intensively  by  Landweber  and  Chwang^10l  Guo  and  Chwang^12^  applied  their  treatment  to  a 
pair  of  two-dimensional  cylinders  and  modified  the  integral  equations  for  accurate  solutions.  For  a 
pair  of  three-dimensional  bodies  translating  in  an  unbounded  fluid,  the  surface  source  distributions 
E(P')  at  point  F  on  surface  1  and  F(Q')  at  point  Q'  on  surface  2  are  governed  by  the  integral  equations 


rtE(P)- Jl 
Si 


E(P')-—  p^-dS 
3Ni  Ril 


2  it  F(Q)  -  fEfPV^-jr-dS 
J  3N2  k12 
Si 


S2 


J-dS  =  U\ 

K21 


— +u2  — . 

3Ni  3Ni 


p^-dS  =  U3 
K22 


— +  u'4-^-, 

9N2  5N2 


(38a) 


(38b) 


where  Rjj  is  the  distance  between  the  source  point  on  body  i  and  the  field  point  on  the  surface  of  body 
j  (ij  =1,2),  Nj  denotes  distance  in  the  outward  normal  direction  at  the  field  point  on  the  surface  of 
body  i.  The  terms  27tE(P)  in  (38a)  and  2tcF(Q)  in  (38b)  appear  after  differentiating  the  singular 
integrals.  As  long  as  two  bodies  are  separated,  equations  (38a,b)  are  well-defined  at  all  points  P  and 
Q.  However,  there  are  two  types  of  difficulties  need  to  be  considered  in  the  numerical  integrations: 

1.  There  are  apparent  singularities  involved  in  the  kernels  of  the  first  integral  in  (38a)  and  of  the 
second  integral  in  (38b)  when  the  point  of  integration  coincides  with  the  field  point  on  the  same 
surface.  As  discussed  by  Landweber  and  Chwang^10^,  the  effect  of  apparent  singularities  can  be 
reduced  significantly  by  noting  that  the  flux  through  a  closed  surface  due  to  a  unit  source  on  the  same 
surface  is  2n.  For  example,  if  there  is  a  unit  source  at  P  on  surface  1,  then 

’  f  R— dS  =  2k  •  (39) 

J  9N  i  Kn 
si 

Thus,  the  first  integral  in  (38a)  can  be  modified  as 
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(40) 


J-dS  =2jiE(P)  - 
Kli 


si 


nr —  E(p> 

Kll 


dN'i 


R 


-]dS, 


11 


where  the  singularity  at  P  =  P  is  eliminated. 

2.  There  exist  steep  peak  values  in  the  kernels  of  the  second  integral  in  (38a)  and  of  the  first 
integral  in  (38b)  when  one  body  is  in  the  proximity  of  the  other  and  the  point  of  integration  is  not  on 
the  same  surface  as  that  of  the  field  point.  Let  8  be  the  gap  distance  between  two  surfaces,  the 
maximum  peak  value  of  a  typical  kernel  in  the  integral  equations  is  of  the  order  of  (1/82)  for  three- 
dimensional  bodies.  In  order  to  remove  the  peaks,  we  may  use  a  similar  treatment  to  subtract  a  term 
from  the  integrand  and  then  to  add  its  accurate  integration  back  to  the  equation.  For  example,  the 
second  integral  in  (38a)  may  be  written  as 


p~dS  = 
K21 


S2 


1F(Q>F(Q0)l^4“'-0' 


^  +  F(Q„) 


S2 


•_a_  1,0 

3N,  ’ 


(41) 


in  which  Qq  is  the  point  where  the  maximum  peak  value  occurs.  The  result  of  the  integration  on  the 
left-hand  side  will  be  improved  if  the  second  integration  on  the  right-hand  side  can  be  obtained  exactly 
or  computed  more  accurately  without  too  much  trouble.  However,  this  modification  has  only 
weakened,  rather  than  removed,  the  ill-behaved  kernel  since  another  peak  value,  whose  magnitude  is 
smaller  than  the  original  one,  is  still  present  in  the  first  derivative  of  the  kernel  on  the  left-hand  side. 
Depending  on  the  smoothness  of  solid  surfaces,  this  modification  may  be  used  successively^17^. 

The  accuracy  of  the  numerical  solution  for  planar  motion  of  two  spheres  will  be  examined  by 
comparing  it  with  the  analytic  solution.  Let’s  define  two  spherical  polar  coordinate  systems  (rt,  a,  (3) 
and  (r2,  0, 8)  by  (see  Fig.  4) 


x’  =  rt  since  cos(3  +  s  =  r2sin0  cos8, 
y'  =  rt  since  sinp  =  r2sin0  sin8, 

z’  =  rj  cosce  =  r2  cos0.  (42) 

From  equations  (38),  (40)  and  (41),  the  surface-source  distribution  on  each  spherical  surface  satisfies 
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the  pair  of  integral  equations 
2k  k 

4jiE(a,P)-a2  Jdp’  {[ECa’.pO-ECa.pJlKuCa.p.a'.pOda' 
o  o 

2k  k 

-b2  Jd8’  J[F(0’,S')-F(O,  £)]K2i(a,p,0',8')d0' 
o  o 

2k  k 

=b2F(0,^)  Jd8'  jK2i(a,p,0',8')d0'  +  U  icospsina+l^sinpsina, 
o  o 
2k  k 

47tF(0,8)-a2  Jdp’  J[E(a’,p')-E(7i,|)]Ki2(9,S,a',p’)da' 
o  o 

2tc  k 

-b2  JdS*  J[F(0',5,)-F(0,8)]K22(0,8,0',5,)de' 
o  o 

2:t  k 

=a2E(7t,^)  JdP'  jKi2(0,5,a',P')da'+U3COs8sin0+U4sin5sin0. 
o  o 


(43a) 


(43b) 


Define 


R($tT1,£  ;  x,y,z)  =  [(^  -  x)2  +  (q  -  y)2  +  (£  -  z)2Y'2 , 


then  the  distance  between  the  source  point  on  body  i  and  the  field  point  on  body  j  is  (Fig.4) 

Ru=  (s+a  sina’cosP’,  a  sina'sinP'  ,a  cosa';  s+a  sinacosP,  a  sinasinp,  a  cosa), 

R21=  (b  sin0'cos8',  b  sin0'sin8',b  cos0';  s+a  sinacosP,  a  sinasinp,  a  cosa), 

Ri2=  (s+a  sina'cosP',  a  sina'sinP',  a  cosa';  b  sin0cos5,  b  sin0sin5,  b  cos0), 

R22=  (b  sin0'cos8\  b  sin0'sin8',  b  cos0';  b  sin0cosS,  b  sin0sin8,  b  cos0),  (44) 


where  s  is  the  separation  distance  between  two  spheres,  s  =  0j02.  The  kernel  functions  in  the  integral 
equations  (43a,b)  are  given  explicitly  by 


Kn(a,p,a',P’) 


_ sina' _ 

23/2  a2  (l-sinasina'cos(p-p')  -  cosacosa')1/2 
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K2i(a,p,0',8')  =  -  -4-  sin0'[a-b  sin0'sinacos(8'-p)-b  cosO'cosa  +s  sinacosP]  , 

Ki2(0,5,a',P')  =  -  -^-sina'[b-a  sina'sin0cos(P'-8)-a  cosa'cos0-s  sinOcosS] , 
Rj2 


K22(e,8,0\8’)  =  - 


23/2  b2  (l-sin0sin0'cos(S-8')  -  cos0cos0')1/2 


In  the  derivation  of  equations  (43a,b),  we  have  used  the  identity  (39)  and  the  relation  for  the  i-th 


sphere 


—  (-L)  =  — (-L) 

3Ni  Rii  3N(  Rii 


Based  on  surface-source  distributions  obtained  from  integral  equations  (43a,b),  we  can 
compute  the  added-mass  coefficients  due  to  the  motion  of  each  sphere  by  means  of  equation  (37)  and 
the  generalized  Taylor's  formula  (21).  For  U'i  =  1,  U'2=U’3=U'4=0, 

3  27t  n 

ku  =  ^3  (3  JdP  J  E(a,P)  sin2a  cosP  da  -  1),  (47a) 

b  o  o 


kn  =  3  Jd8  J  F(0,8)  sin20  cos8  d0. 


(47b) 


o  o 


For  U'2  =  1,  U'1=U,3=U’4=0, 


3  271  K 

k22  =  ^j(3  Jd  P  J  E(a,P)  sin2a  sinP  da  -  1), 
b  o  o 


k24  =  3  Jd 8  J  F(0,8)  sin20  sin8  d0. 
o  o 

Similarly,  for  U’3=l,  U'i=U,2=U’4=0, 


k33  =  3  Jd8  J  F(0,8)  sin20  cosS  d0  -1. 
o  o 


(47d) 
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Finally,  for  U’4=l,  U'i=U’2=U’3=0, 


IC44  =  3  JdS  J  F(0,S)  sin20  sin8  d0  -  1. 
o  o 

As  explained  in  deriving  equation  (1 1),  we  have 
^12  =  ^23  =  ^14  =  ^34  =0- 


(470 


(47g) 


IV.2  Numerical  solutions  of  integral  equations 

The  unknown  source  distributions  E(a,P)  and  F(0,5)  for  two  spheres  are  governed  by  a  set  of 
two  integral  equations  (43a,b).  These  equations  will  be  transformed  into  two  sets  of  linear  equations 
by  introducing  appropriate  quadrature  formulas  for  the  integrals,  and  solved  by  the  Gauss-Seidel 
iterative  method.  Thus,  for  the  (n+l)-th  iteration,  we  have 


4*E"*‘  -  a2S  £  W.;W,,|F.',  -  EJ]  KlltjiT  -  b2 1  I  Qlk.Q2n,[F;.m.  -I*]  K21ijk.m 


1  J 

L 

k’  m 


k’  m’ 


=  b2F^XI  i2lk.Q2m.  K21ijk.m.  +  U1!  cosp  sina  +  U'2  sinP  sina. 


(48a) 


‘‘"C-a2?  I  K12tmirb2 17  r  lkfi2m  K22kmk  m. 

i  j  4  4  1  m 


2n+l  r,n+l 


=  a2  Eg+1  X  X  wii,w2j'  K12kmi  j’  +  ^3  cos$  sin0  +  U'4sin8  sin0. 


(48b) 


where  W,  £2  are  weighting  factors  corresponding  to  different  quadrature  formulas,  subscripts  i,  j,  k,  m 
refer  to  the  location  of  a  field  point  along  a,  p,  0,  5  directions  respectively,  i',  j',  k',  m'  refer  to  that 
of  a  source  point,  and  the  superscript  n  stands  for  the  n-th  iteration.  To  start  the  numerical  iteration 
process,  we  assume  that  the  first  approximation  is 

=  (U'j  cosP  sina  +  U'2sinP  sina)/47t  (49a) 


and 


Fkm  =  (U’3  cos8  sin0  +  U’4sin8  sin0)/47t.  (49b) 

In  fact,  based  on  the  numerical  results,  equations  (48a)  and  (48b)  always  converge  to  the  correct 
solution  regardless  of  the  initial  approximation  (49). 
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When  two  spheres  are  close  to  each  other,  it  is  evident  that  the  source  distribution  on  each 
surface  changes  rapidly  around  the  gap  between  two  surfaces.  A  typical  distribution  of  E(a,p)  for  the 
case  where  ITj  =  1,  U'2  =  U’3  =  U'4  =  0,  and  a/b  =  1.0,  s/b  =  2.03  is  shown  in  Figure  5  .  To  obtain 


an  accurate  numerical  integration  over  the  spherical  surface,  we  shall  put  more  node  points  near  the  gap 
by  subdividing  the  region  of  integration  into  few  partitions  and  apply  the  Gaussian  quadrature  formula 
in  each  subregion.  The  sizes  of  partitions  measured  by  coa,  cop,  etc.  (Fig.6)  are  dependent  on  the 

separation  distance  s  of  two  spheres.  It  is  well  known  that  most  quadrature  formulas,  including  the 
Gaussian  quadrature  formula,  performs  the  best  when  the  ill-behaved  point  is  at  the  boundary  of  the 
integration  region.  In  solving  the  integral  equations  (48a,b),  however,  we  need  to  fix  the  field  points 
as  well  as  the  point  of  integration  on  each  surface  in  order  to  iterate  the  unknown  source  distributions. 
Thus,  to  compute  the  integrations  on  the  right-hand  sides  of  (43a,b)  accurately,  we  shall  rotate  the 
coordinates  in  such  a  way  that  for  each  fixed  field  point  on  one  surface,  the  maximum  peak  value  on 
the  other  is  always  at  the  boundary  of  the  integration  region.  Suppose  that  the  peak  of  K2  j , 
corresponding  to  a  field  point  P(a,  a0,  P0)  on  sphere  1,  occurs  at  q(b,  0*,  5*)  on  sphere  2  (Fig.  7), 
then  the  integration  of  K21  over  surface  2  should  be  performed  in  the  (X',  Y',  Z')  coordinate  system 
which  is  obtained  by  first  rotating  the  (x\  y’,  z')  coordinate  system  about  the  z'  axis  with  an  angle  5*. 
then  rotating  about  the  new  y'  axis  with  an  angle  (j  -  0*).  Therefore, 


T  X’, 


where 


T  = 


cos5*sin0* 
sin5  sin0 
COS0* 


-sin5*  -cos8*cos0* 
cos8*  -sin5*cos0* 
o  sin0* 


(50a) 


(50b) 


With  this  coordinate  transformation,  the  first  and  the  last  abscissas  for  the  Gaussian  quadrature 
formulas  in  both  0  and  5  directions  are  in  the  nearest  neighbourhood  of  q. 


However,  the  explicit  expression  for  point  q  in  terms  of  a0  and  (30  is  complicated  so  that  we 
may  have  recourse  to  an  approximation.  It  is  noted  that  the  peak  affects  the  numerical  integration 
significantly  only  in  a  region  when  the  gap  between  two  spheres  are  very  small  and  the  field  point  is 


around  the  centerline  (Fig.  7).  In  this  region,  the  location  of  q^  which  is  on  the  line  of  Po2  (Fig.  7), 
should  be  very  close  to  q,  where  the  peak  value  of  K21  takes  place.  Thus 


T,  ,sin8\  „  „  (1  -  y2/2)  sin0' 

K21=  -(-J-)  COS(7t-V) - -2 - 

*21  **21 


(V«  1). 


We  may  examine  the  above-mentioned  improvement  by  considering  the  following  integral  over  sphere 
2: 


2k 


T  fJ5,t  f  sin6’  Jrtl  in  ,  ,d-b, 

I=  JdS  (  -^-de  =-rdln(dTb)- 

o  o  *21 


2ti  ,  ,d-b, 


(51) 


where  d  =  IPo2l  (Fig.7).  In  this  integration,  the  position  of  the  source  point  q  is  exactly  on  the  line  Po2. 
The  analytical  and  numerical  results  of  the  integration  I  are  plotted  in  Fig.  8,  where  we  note  that  the 
deviation  of  the  modified  numerical  result  from  the  exact  one  is  much  smaller  than  that  corresponding 
to  the  unmodified  numerical  result.  Applying  the  modification  (50)  to  equation  (48a,b)  and  by  means 
of  (47),  we  can  compute  the  added-mass  coefficients.  When  a/b  =  1.0  and  s  varies  from  2.02  to  10, 
the  numerical  results  of  added-mass  coefficients  are  listed  in  Table  2. 


V.  DISCUSSION  OF  NUMERICAL  RESULTS 
The  Lagrange's  equations  of  motion,  associated  with  the  expression  of  kinetic  energy  of  the 
fluid  due  to  the  planar  translation  of  two  bodies  of  revolution,  can  be  applied  to  determine  the 
hydrodynamic  interaction  between  two  bodies,  as  long  as  the  flow  is  a  potential  flow.  To  solve 
equations  (10)  explicitly,  we  shall  start  the  computation  at  the  initial  position  xjo  (i=l,  2,  3, 4)  with  the 
initial  absolute  velocity  components  uio  and  the  external  forces  Ej.  At  the  beginning  of  the  j-th  time 
interval  8tj  (j  =  1,  2,  ...),  the  absolute  velocities  Ujj  are  obtained  by  the  fifth-order  Runge-Kutta 
integration  in  terms  of  velocities  in  the  previous  time  interval,  u^.j,  and  the  added  masses  and  their 
derivatives  evaluated  at  Xjj.  At  the  end  of  the  j-th  time  interval,  the  new  position  of  each  body  is 
simply  calculated  by  Xjj  =  Xjj.,  +  u^Stj.  The  size  of  5tj  is  adaptive  according  to  the  value  of  a  pre¬ 
assigned  error-control  parameter  (1.0*10’5)  and  the  difference  between  the  fifth-  and  the  sixth-order 
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Runge-Kutta  integration  of  u^.  This  process  is  repeated  until  either  two  bodies  are  in  contact  or  one 
body  passes  over  the  other. 

There  are  a  number  of  real  geometrical  situations  in  which  one  can  apply  the  afore-mentioned 
mathematical  model  to  predict  the  motion  of  solids.  In  the  present  study,  both  bodies  are  assumed  to 
be  spheres  with  radii  a  and  b  (a<b)  for  bodies  1  and  2  respectively.  In  addition,  we  shall  normalize  all 
lengths  by  b,  a/b  <,  1.0  and  s/b  >  (1.0+a/b). 

We  first  consider  the  motion  of  a  spherical  particle  of  radius  a,  conveyed  by  a  uniform  flow, 
around  a  large  spherical  body  fixed  in  space.  This  problem  has  applications  in  the  ice-coating  process. 
With  a/b=0.1,  xlo/b  =  -20,  x2o/b  =  0.1,  0.3,  etc.,  and  ul0/U0  =  1.0,  u2o/U0  =  0.0,  where  U0  is  the 

uniform  flow  in  the  x  direction,  the  trajectories  of  the  particle  are  plotted  in  Fig.  9a  to  9e  for  different 
density  ratios  of  the  body  to  the  fluid  medium.  Fig.  9a  corresponds  to  the  case  of  an  ice  particle 
moving  in  fresh  water  and  Fig.  9e  shows  the  same  ice  particle  in  an  air  flow.  The  ratios  2.0,  5.0,  and 
10.0  correspond  to  various  fluid  media  and  illustrate  the  change  of  trajectories  with  respect  to  the 
density  ratio.  From  these  figures,  we  observe  that  the  motion  of  the  particle  is  affected  by  its  inertia, 
which  prevents  the  particle  from  moving  out  of  its  straight  path,  and  by  the  interaction  with  the  second 
sphere,  which  bends  the  particle  trajectory  to  a  curved  streamline.  In  the  case  of  an  ice  particle  carried 
by  an  air  flow  (Fig.  9e),  the  inertia  effect  is  so  predominant  that  the  trajectories  are  almost  straight 
lines,  whereas  for  the  same  particle  moving  in  water  (Fig.  9a),  the  curvature  of  the  trajectories 
becomes  very  large  when  two  bodies  are  close  to  each  other.  This  conclusion  can  also  be  drawn  from 
the  equations  of  motion  (20),  which  indicate  that  the  accelerations  of  the  particle  in  both  x,  and  x2 
directions  are  proportional  to  (p^p)'1.  In  some  physical  problems,  it  is  important  to  determine 
whether  or  not  a  drifting  body,  conveyed  by  a  current,  can  impact  with  a  fixed  body.  This  physical 
property  can  be  expressed  by  a  "collection  coefficient"  E  which  is  defined  as  the  ratio  of  the  critical 
initial  position  x*^,  below  which  the  drifting  body  will  impact  with  the  fixed  body,  to  the  radius  of  the 
fixed  body.  Fig.  10  shows  the  result  of  the  collection  coefficient  for  a  pair  of  spheres,  0.01  <  a/b  <  1 . 
in  fluids  of  different  densities.  Moving  bodies  of  very  small  size  are  not  considered  since  the 
Reynolds  number  becomes  quite  small  such  that  the  inviscid-fluid  theory  is  not  applicable. 
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Regarding  a  floating  body  in  sea  water,  on  the  other  hand,  we  should  consider  the  change  of 
hydrodynamic  interaction  between  an  ice  floe,  idealized  as  a  sphere,  and  a  fixed  spherical  offshore- 
structure  versus  the  size  of  the  floating  body.  For  p^p  =  0.89,  the  trajectories  of  a  floating  sphere  of 

radius  a  (a/b=0.5,  1.0)  around  a  fixed  sphere  are  shown  in  Fig.  11.  Initially,  the  floating  sphere 
moves  with  the  uniform  flow,  ul0/U0  =  1.0  and  u2o/U0  =  0.0.  These  trajectories  are  quite  flat  in 

comparison  with  the  ones  of  a  small  body  since  the  inertia  effect  of  a  large  body  predominates  over  the 
hydrodynamic  interaction  force  due  to  the  presence  of  a  second  body. 

To  illustrate  the  dependence  of  trajectories  of  a  moving  sphere  on  its  size,  we  plot  the  velocity 
components  UjAJ0  and  u2/U0  in  Fig.  12,  and  the  trajectories  in  Fig.  13,  respectively,  for  a  sphere  of 

various  radius  ratios  a/b  from  0.1  to  1.0  around  a  fixed  sphere  of  radius  b.  The  initial  conditions  used 
in  these  two  figures  are  xlo/b  =  -20,  x2o/b  =  0.5,  ulo/Uo=1.0,  U2cA^o=  0*0*  a  f"lxed  density  ratio 
p,j/p=0.89.  When  two  bodies  are  close  to  each  other.  Fig.  12  shows  slight  differences  on  the  velocity 
components  due  to  the  size  variation,  while  Fig.  13  shows  that  for  the  same  initial  position,  the 
trajectories  of  sphere  1  is  almost  independent  of  its  size.  In  Fig.  13,  only  the  trajectory  for  the  sphere 
a/b=0.1  (dotted  line)  can  reach  the  position  xt  =  0,  which  indicates  that  the  small  sphere  passes  over 
the  fixed  one.  All  other  trajectories  terminate  at  certain  position  o at  which  the  two  spheres  are  in 

contact.  This  phenomenon  can  be  understood  by  considering  the  added  masses.  The  first  terms  of 
added  masses  in  (22)  and  (34)  and  their  derivatives  in  (23)  and  (35)  are  proportional  to  a3.  These 
terms  are  predominant  when  s  is  sufficiently  large,  while  the  rest  becomes  significant  only  when  s  is 
very  close  to  (a+b).  Consequently,  the  predominant  parts  of  the  solution  of  equation  (10)  are  almost 
independent  of  the  radius  a,  since  a3  appears  on  both  sides  of  (10).  Based  on  this  observation,  we  can 
surely  use  the  trajectories  for  a  small  particle  around  a  fixed  blunt  body  as  a  good  first  approximation 
for  the  trajectories  of  a  large  body,  or  in  some  cases,  we  can  even  neglect  the  effect  due  to  the  size 
variation. 

In  the  case  of  a  small  body  pursuing  a  large  target  which  moves  in  a  fluid,  we  shall  consider 
the  influence  of  the  motion  of  the  target  on  the  trajectories  of  the  small  body.  Suppose  that  a  large 
spherical  target  of  radius  b  with  pjp  =  0.89  is  initially  located  at  (10,0)  and  moves  in  a  stationary  fluid 
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with  u3o=1.0,  u^sO.O.  A  small  panicle  of  radius  a  (a/b=0.1)  and  density  pJp=Q&9  is  released  at  (0, 
x^)  with  initial  velocities  ul0  =  2.0  and  u2o  =  (x^-l.lj/lO.  We  can  determine  its  trajectories  based  on 
the  solution  of  (20).  Without  interaction,  we  would  expect  that  these  two  bodies  come  in  contact  at 
x=20.  However,  from  Fig.  14,  we  note  that  the  trajectories  of  the  small  panicle  are  bent  significantly 
in  a  region  close  to  the  large  target  due  to  the  hydrodynamic  interaction.  This  plot  is  consistent  with 
the  physical  interpretation. 

The  magnitude  and  the  direction  of  the  hydrodynamic  interaction  force  between  two  bodies 
depend  on  the  relative  motion  between  them  and  the  direction  of  the  oncoming  flow.  Fig.  15  shows 
the  variation  of  the  magnitude  of  the  interaction  force  versus  the  flow  direction  and  the  separation 
distance.  We  note  from  Fig.  15  that  if  the  angle  between  the  centerline  and  the  flow  direction  is  less 
than  a  certain  value  depending  on  the  separation  distance,  the  interaction  force  is  repulsive,  which 
prevents  the  collision  of  two  bodies,  whereas  if  the  angle  is  larger  than  this  value,  the  interaction  force 
becomes  attractive.  Thus,  when  a  particle  moves  towards  the  target  sidewise  from  behind  (see  Fig. 
14)  and  the  angle  between  the  centerline  and  the  oncoming  current  is  small,  the  interaction  force  repels 
the  particle  from  the  target  and  reduces  its  velocity  component  u2.  On  the  other  hand,  when  the  particle 

moves  alongside  of  the  target,  the  relative  oncoming  flow  becomes  almost  perpendicular  to  the 
centerline,  and  the  interaction  force  pushes  the  particle  toward  the  target.  This  conclusion  also 
indicates  that  if  the  angle  between  the  centerline  and  the  relative  flow  is  large,  the  motion  of  the  target 
generates  an  attractive  force  and  the  trajectory  of  the  particle  points  toward  the  target  without  much 
bending  (Fig.  14,  x2o/b  =  4.0). 

Another  practical  problem  we  shall  consider  in  the  present  study  is  the  hydrodynamic 
interaction  between  two  bodies,  idealized  as  spheres,  when  they  move  arbitrarily  in  a  stationary  fluid. 
Fig.  16  shows  trajectories  of  two  equal  spheres  (a/b  =  1)  with  pjp  =  0.89,  which  corresponds  to  ice 
floes  in  sea  water.  Initially,  both  sphere  move  in  a  stationary  fluid  with  velocities  Uj  =  u3  =  1,  u2  =  u4 
=  0,  and  the  initial  separation  distance  sjb  varies  from  2.2  to  3.0.  As  we  have  discussed  previously, 

the  interaction  force  due  to  this  type  of  motion  attracts  both  bodies  and  drives  them  toward  each  other. 
The  strength  of  the  force,  shown  in  Fig.  17  with  a/b  =  1  and  u,  =  u3  =  1,  surely  depend  on  the 
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separation  distance  s/b.  The  time  required  for  two  equal  spheres  to  come  into  contact  is  shown  in  Fig. 
18  versus  the  initial  separation  distance  Sq/I).  For  pjp  =  0.89  and  a/b  =  0.3,  0.5,  0.7,  1.0,  Fig.  19 
shows  the  trajectories  of  two  spheres  when  they  are  initially  released  at  xl0  =  x3o  =  0  and  x2o  =  -x^  = 
1.25b  with  velocities  ulo  =  u3o  =  1  and  u2o  =  u^  =  0  in  a  stationary  fluid.  From  this  figure,  we  note 
that  if  sphere  1  is  much  smaller  than  sphere  2,  its  motion  cannot  affect  the  motion  of  sphere  2 
significantly,  since  the  inertia  of  sphere  2  is  much  larger  than  that  of  sphere  1.  We  also  note  that  the 
variation  of  trajectories  of  sphere  1  due  to  the  variation  of  the  radius  ratio  a/b  is  much  smaller  than  that 
for  sphere  2,  since  the  added  masses  associated  with  sphere  2  are  more  sensitive  to  the  change  of  a/b. 

VI.  CONCLUSIONS 

A  potential-flow  prediction  for  the  motion  of  a  pair  of  bodies,  and  the  hydrodynamic  interaction 
force  acting  upon  them  in  an  inviscid  fluid,  have  been  presented.  The  Lagrange's  equations  of  motion 
are  generalized  for  planar  translational  motions  of  bodies,  including  the  effects  of  solid  constraints, 
external  forces  in  the  plane  of  motion,  and  a  uniform  stream  in  any  direction  parallel  to  the  plane  of 
motion.  This  generalization  is  applicable  to  bodies  of  revolution  which  are  symmetric  with  respect  to, 
and  have  their  axes  of  rotation  perpendicular  to,  the  plane  of  motion.  The  velocity  components  and  the 
moving  trajectories  of  each  body  are  obtained  by  integrating  the  equations  of  motion  in  terms  of  the 
kinetic  energy  of  the  fluid  which  is  a  function  of  added  masses. 

In  order  to  determine  the  reliability  of  numerical  solutions  of  interactions  between  a  pair  of 
three-dimensional  solids,  the  exact  solution  of  added  masses  and  their  derivatives  in  closed  forms  is 
first  considered  for  the  centroidal  and  transversal  motions  of  two  spheres.  A  new  iterative  formula  has 
been  developed  to  evaluate  the  added  masses  due  to  the  transversal  motion  of  two  spheres  based  on  the 
analysis  of  the  velocity  potential  in  the  near  field  around  the  bodies.  The  added  masses  due  to  the 
centroidal  motion  are  obtained  by  means  of  the  Taylor's  added-mass  formula  and  an  iterative  scheme 
for  the  strengths  and  positions  of  interior  doublets. 

The  boundary- integral  method  and  the  generalized  Taylor's  added-mass  formula  are  used  for 
numerical  solutions  of  added  masses.  However,  when  two  bodies  are  very  close  to  each  other,  it  is 
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difficult  to  obtain  an  accurate  result  because  of  the  ill-behaved  integral  equations.  These  integral 
equations  are  modified  by  first  subtracting  out  the  steep  peak  value  from  the  integration  so  that  the 
resultant  kernel  is  bounded,  and  then  adding  an  accurate  integration  of  the  kernel  back  to  the  equation. 
As  an  example,  the  added  masses  of  two  equal  spheres  are  computed  numerically  and  the  results  are 
compared  with  exact  solutions.  A  very  good  agreement  has  been  obtained. 

Due  to  the  limitation  of  the  potential-flow  analysis,  the  afore-mentioned  results  are  only 
applicable  to  cases  where  the  effects  of  fluid  inertia  and  the  nonuniformity  of  the  flow  due  to  the 
presence  of  a  second  body  are  dominant. 
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Table  1.  Exact  Values  of  Added-Mass  Coefficients 
for  Two  Equal  Spheres  fa/b  =  LO) 


^11 

^13 

.^22 

dk24 

s 

kll 

ds 

"  k13 

ds 

k22 

ds 

*  k24 

ds 

2.01 

0.5702 

0.5708 

0.2163 

0.7257 

0.5191 

0.1149 

0.0984 

0.2078 

2.02 

0.5651 

0.4486 

0.2097 

0.6020 

0.5180 

0.1005 

0.0964 

0.1923 

2.03 

0.5610 

0.3796 

0.2041 

0.5315 

0.5171 

0.0901 

0.0945 

0.1808 

2.04 

0.5575 

0.3321 

0.1990 

0.4825 

0.5162 

0.0819 

0.0927 

0.1716 

2.05 

0.5543 

0.2963 

0.1944 

0.4453 

0.5154 

0.0752 

0.0911 

0.1638 

2.06 

0.5515 

0.2679 

0.1901 

0.4154 

0.5147 

0.0695 

0.0895 

0.1571 

2.07 

0.5490 

0.2445 

0.1861 

0.3906 

0.5140 

0.0646 

0.0879 

0.1511 

2.08 

0.5466 

0.2248 

0.1823 

0.3695 

0.5134 

0.0603 

0.0864 

0.1458 

2.10 

0.5425 

0.1930 

0.1752 

0.3349 

0.5123 

0.0530 

0.0836 

0.1365 

10.0 

0.5000 

0.0000 

0.0015 

0.0005 

0.5000 

0.0000 

0.0008 

0.0002 

Table  2.  Numerical  Results  of  Added-Mass  Coefficients 

for  Two  Eaual  Soheres  fa/b 

s 

kn 

'  ki3 

k22 

k24 

“a 

2.02 

0.5614 

0.2063 

0.5182 

0.0967  20 

15 

2.03 

0.5591 

0.2025 

0.5171 

0.0947  20 

15 

2.04 

0.5562 

0.1981 

0.5162 

0.0929  20 

15 

2.05 

0.5540 

0.1943 

0.5154 

0.0910  25 

20 

2.06 

0.5512 

0.1900 

0.5147 

0.0895  25 

20 

2.07 

0.5486 

0.1860 

0.5140 

0.0879  30 

25 

2.08 

0.5464 

0.1823 

0.1534 

0.0864  30 

25 

2.10 

0.5421 

0.1752 

0.5123 

0.0836  30 

35 

10.0 

0.5000 

0.0015 

0.5000 

0.0008  60 
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Figure  Captions 

Figure  1.  Relative  rectangular  coordinate  system  (x,  y,  z)  and  velocity  components. 

Figure  2.  Relative  rectangular  coordinate  system  (x1,  y',  z')  and  velocity  components. 

Figure  3.  Spherical  polar  coordinate  systems  (rlt  0j,  Xj)  and  (r2,  02,  X2). 

Figure  4.  Surface  integration  of  source  distributions  on  spheres  1  and  2. 

Figure  5.  Three-dimensional  plotting  of  surface  source  distribution  on  sphere  1  with  a/b  =  1 ,  s/b  = 
2.03,  U\=l,  and  U2=  U3=  U4=  0. 

Figure  6.  Definition  of  regions  containing  peaks. 

Figure  7.  Rotation  of  the  coordinate  system  (x\  y',  z')  and  definition  sketch  for  d,  a0,  and  (3Q. 

Figure  8.  Comparison  of  numerical  result  of  integral  I  (51)  obtained  with  and  without  rotating 
coordinates. 

Figure  9.  Trajectories  of  a  moving  sphere  in  different  fluids  with  a/b=0.1,  ulo/U0=l,  u2o/U0=  0, 
and  xlo/b  =  -20. 

Figure  10.  Dependence  of  the  collection  coefficient  E  of  an  ice  particle  on  the  size  for  various  density 
ratios. 

Figure  11.  Trajectories  of  .  *  .oving  sphere  with  xl0/b  =  -20,  ulo/U0=l,  u2o/U0=  0,  and  pa/p  = 
0.89. 

Figure  12.  Velocity  components  ut  and  u2  of  a  moving  sphere  with  various  radius  ratios. 

Figure  13.  Trajectories  of  a  moving  sphere  with  various  radius  ratios. 

Figure  14.  Trajectories  of  a  moving  sphere  with  radius  a  affected  by  another  sphere  with  radius  b, 
a/b=0.1,  p^p  =  0.89. 

Figure  15.  Dependence  of  interaction  forces  on  the  direction  of  a  uniform  flow. 

Figure  16.  Trajectories  of  two  equal  spheres  (a/b=l)  moving  initially  perpendicular  to  their  centerline 
in  a  stationary  fluid  with  pa/p  =  0.89,  a/b=l,  xl0/b  =  x3o/b  =  0,  ulo  =u3o  =  1,  and  u-,0 

=u4o  =  0. 

Figure  17.  Interaction  force  between  two  equal  spheres  (a/b=l )  versus  the  separation  distance  s/b. 

Figure  18.  Time  required  for  two  equal  spheres  (a/b=l)  to  come  into  contact  versus  the  separation 
distance  s^b  with  p^p  =  0.89,  xlo/b  =  x3o/b  =  0,  ulo=u3o=  1,  and  u2o=u4o=  0. 

Figure  19.  Trajectories  of  two  spheres  moving  initially  parallel  to  each  other  with  pa/p  =  0.89,  x]o/b 
=  x3o/b  =  0,  x2o/b  =  -  xjb  =  1 .25,  u  1o=u3o=  1 ,  and  u^u^  0. 


Figure  1.  Relative  rectangular  coordinate  system  (x,  y,  z)  and  velocity  components. 


Figure  2.  Relative  rectangular  coordinate  system  (x',  y’,  z')  and  velocity  components. 
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Figure  17.  Interaction  force  between  two  equal  spheres  (a/b=l)  versus  the  separation  distance  s/b. 
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